{-# OPTIONS -XTypeOperators #-}

module Instrument.PSF where

import HSimMath ( (*>), (<->) )
import FTable   ( (:->:)(..), fromList, fromList2, liftFT, liftFT2, extrapolate, move2, scale2 )

import qualified Data.Array as A
import qualified ZEMAX      as Z


type Value = Int :->: (Int :->: Double)

-- точечный источник света

data PointSource = PointSource SED (Double,Double)

type SED = Int :->: Double

-- Представление ФРС (функции рассеяния точки).

type PSF = (Double, Double) -> Int :->: Value

-- Моделирование дисперсии 

disperse :: Double -> Double
-- disperse l = 300  * (l - 1000) / (6000-1000)
disperse l = 0.06 * (l - 1000)

distort :: (Double,Double) -> (Double,Double)
distort (x,y) = (x',y')
    where
       x' = (x + 210) * k
       y' = (y + 210) * k
       k  = 1096/410


-- Создает ФСР в виде 2-ух мерной функции гаусса

mkGaussPsf :: Double -> Double -> PSF
mkGaussPsf sx sy = psf
   where
      psf (u,v) = liftFT (1000,6000,5) (psf' sx sy (u,v))

psf' :: Double -> Double -> (Double, Double) -> Int -> Value
psf' sx sy (u,v) l = liftFT2 (a,b) (\x y -> gauss2d (fromIntegral x + 0.5 - u-d, fromIntegral y + 0.5 -v) )
   where
      a = (floor   (d+u-rx), floor   (v-ry), 1)
      b = (ceiling (d+u+rx), ceiling (v+ry), 1)

      d = disperse $ fromIntegral l

      (rx, ry) = (sx, sy) *> 5

      gauss2d (x,y) = exp( -0.5*( (x/sx)^2 + (y/sy)^2) ) / ( 2 * pi * sx * sy )


-- Загрузка ФРС из файлов, создаваемых программой ZEMAX 

loadPSF :: (Int, Double) -> [FilePath] -> IO PSF
loadPSF i fs = do
  d <- (sequence . map Z.loadPSF) fs
  case sequence d of
    Left  e -> return $ error e
    Right d -> return $ genPSF d i


genPSF ds (h,s) = psf
  where
     -- ФРТ зависит от положения (x,y) точечного источника относительно оси телескопа.
     -- Поэтому ФРТ - это функция действительных значений (x,y).
     psf (x,y)  = liftFT (1000,6000,5) (psf' x y)

     psf' x y l = scale2 (s,s) $ move2 (dx,dy) $ (extrapolate.extrapolate ft) x y
        where
           (dx, dy) = (x', y') *> (1/s) <-> (17,17)
           (x', y') = distort (x,y)
           d        = disperse $ fromIntegral l
           
           ft = fromList2 ((-h,h,h), (-h,h,h)) $ map prepare ds
              where
                 prepare d = ((round $ x*1e4, round $ y*1e4), ft/19)
                    where
                      (x,y)= Z.centerDeg d
                      arr  = Z.psfArray  d

                      ((a,b), (a',b')) = A.bounds arr

                      ft = liftFT2 ((a,a',1),(b,b',1)) (\x y -> arr A.! (x,y))